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Abstract 

A spatial three level food web model with a closed nutrient cycle is presented 
and analyzed via Monte Carlo simulations. The time evolution of the model 
reveals two asymptotic states: an absorbing one with all species being extinct, 
and a coexisting one, in which concentrations of all species are non-zero. 
There are two possible ways for the system to reach the absorbing state. In 
some cases the densities increase very quickly at the beginning of a simulation 
and then decline slowly and almost monotonically. In others, well pronounced 
peaks in the R, C and D densities appear regularly before the extinction. 
Those peaks correspond to density outbursts (waves) traveling through the 
system. We investigate the mechanisms leading to the waves. In particular, 
we show that the percolation of the detritus (i. e. the accumulation of 
nutrients) is necessary for the emergence of the waves. Moreover, our results 
corroborate the hypothesis that top-level predators play an essential role in 
maintaining the stability of a food web (top-down control). 
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1. Introduction 

Every year significant efforts are directed to conservation and restoration 
of natural ecosystems. In order to develop efficient protection strategies of 
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those systems identification and understanding of meclianisms tliat lend tliem 
tlieir stability are required. 

Starting with the pioneering work of Elton f]\, food webs are a central 
organizing theme in studying ecosystems. They depict feeding relationships 
among species in ecological ecosystems. Since their introduction by Elton, 
they have been subject of intensive studies from both theoretical j2|, y, 0, [sl, 0] 
and experimental jl, 0] perspective. 

Although recent decades have seen a significant progress in understand- 
ing of food web properties, our knowledge remains sketchy. Until 1970s the 
predominant belief of ecologists was that large and highly complex ecosys- 
tems were more stable that the simpler ones [if, Isj. It was then questioned by 
May jof, who showed mathematically that large and complex random com- 
munities are inherently unstable. May's have initiated a debate on stability 
of ecosystems which lasts to this day and is still far from being over. 

Several mechanisms have been identified as factors leading to stability of 
food webs. One of them is compartmentalization, i.e. existence of subsets 
of species in the food web that interact more frequently among themselves 
than with other species in the system. It was shown 10| that compartments 
are benefitial to an ecosystem, because they act to buffer the propagation of 
extinctions. A special case of compartments are groups of consumers associ- 
ated with each particular plant species, called component communities 11 . 
A food web consisting of such communities is stable, because a disturbance 
associated with fluctuactions in species' density is confined largely to that 
species' component community. 



Generalist consumers constitute another stabilizing factor [12|. A gen- 
eralist is able to switch from one food source to another, which is more 
abundant. This switching tends to keep a food web stable, since it allows to 



control the abundant species and lets the less common one recover [13 



A process that often occurs in food webs is a top-down control of lower 



trophic levels by apex predators [lj|. The top-level predators play a key 



role in controlling the population of prey and, as a consequence, they limit 
the degree to which the prey endanger primary producers. Many population 
collapses may be traced down to altered top-down forcing regimes associated 
with the loss of native apex predators or the introduction of exotic species. 

Detritus consisting mostly of bodies or fragments of dead organisms has 
been long recognized as one of the important factors in ecology. However, the 
theories of food webs and trophic dynamics have largely neglected detritus- 



based chains and have focused almost exclusively on grazing food chains [15 
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The "green-world" view of these works has been severely criticized 16| and it 
is expected that only merging this approach with the one considering also the 



detritus could yiels a satisfactory ecological theory [15[. The existing theoret 



ical studies [l7| indicate that the flow of energy from detritus to living chain 



may increase the extinction cascades. There are also arguments [18| that 
consumption of prey from detritus chain may weaken top-down regulations 
of stability. Thus, many fundamental questions asked in ecology, concerning 
the structure of food webs, the length of food chains, the size and direction 
of extinction cascades, may have a different form and interpretation when 
detritus is taken into account. 

Food webs consist of species using dispersal as a mean of both procuring 
key resources and avoiding natural enemies. Therefore they are inherently 
spatial entities, a fact which is also neglected in the vast majority of studies 



leading to a rather poor understanding of spatial food web dynamics [19 
An ubiquitous feature of different spatial systems in nature are waves travel- 
ing through them. Those waves are of big importance, because they control 



the speed of many dynamical processes including chemical reactions [20 



epidemic outbreaks [2l| and biological evolution [22]. Although great effort 



has been expended into the description of the waves [22|, |23| , their rigorous 



description is still lacking mainly due to their sensitivity to different kinds 
of fluctuations. Thus, any contribution that provides insight into the mech- 
anisms governing the traveling waves is valuable to our understanding of 
complex systems with spatial structures. 



In our recent study 



2J, |25| a three species food web model with a detri- 



tus path was analysed by means of Monte Carlo simulations. Our findings 
indicate that under certain conditions complex spatio-temporal patterns in 
form of density waves appear in the food web. Analysis of those waves is the 
main goal of this work. 

The paper is organized as follows. In Section [2] the model is introduced. 
Simulation results are discussed in Section [31 And finally in Section H] con- 
clusions are drawn. 



2. Model 

Fig. [T] illustrates a hypothetical system we are going to investigate. The 
food web consists of three trophic levels. The basal level species will be 
called resources (R) within our model. It corresponds to primary producers 
or autotrophs in real ecosystems, i.e. to organisms able to take energy from 
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Figure 1: Model food web investigated in this paper. The arrows indicate the flow of 
nutrient in the trophic loop. Consumers C feed on the resources R and are themselves 
food for predators P. Detritus D consisting of dead consumers and/or predators provides 
nutrient for resources. 



the environment (sunlight or inorganic chemicals) and convert it into energy- 
rich molecules such as carbohydrates. 

The species at intermediate level in Fig. [1] will be called consumers (C). 
It relates to herbivores in real systems, which principally feed on primary 
producers. 

The consumers themselves constitute food for the top level species called 
predators (P). In general, the predators correspond to carnivores in real 
systems. 

The remains of the consumers and predators form detritus (D), which 
provides nutrient for the resources. 

The model food web was built up under the tacit assumption that the 
conversion of dead fragments into nutrient occurs immediately and without 
any external help. In reality, dead organisms are broke down and converted 
into useful chemical products by decomposers. For the sake of simplicity we 
will neglect them in the present work and address their role for the viability 
of the whole food web in a forthcoming paper. 

It should be noted that the model shown in Fig.[T]resembles some real ma- 
rine food webs with phytoplankton being the primary producer, zoo plankton 
and fish occu pyin g the higher trophic levels and bacteria playing the role of 



decomposers [26 



We are going to investigate the model by making use of an agent-based 
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Monte Carlo simulation. To this end we put individuals of each species on a 
square lattice. Each individual is characterized by two parameters: a death 
rate di and a birth rate hi. The death rate determines the probability that 
an organisms dies in a given time step if it is not able to feed. The birth rate 
stands for the ability of an individual to convert food reserves into breeding 
success. Theoretically, at each trophic level there may live different species. 
We will assume however that each level is occupied by only one species. 
Moreover, to keep the model as simple as possible, we do not differentiate 
between organisms of a given species. In other words, each individual of the 
species i {i = R,C or P) will have exactly the same pair of the parameters 
{bi,di). 

Once put randomly on the lattice, the system will evolve due to the 



following rules [24, |25 



While populating the nodes we are applying an exclusion principle stat- 
ing that each lattice node may be occupied by at most one agent of 
a given type. Multiple occupancy of nodes by individuals of different 
species is allowed. 

If a C appears at a node populated by an R, then R is eaten and 
removed from the lattice. This may correspond to a herbivore eating 
a plant. Next, a progeny of C is created with a probability be at a 
randomly chosen node in its Moore neighborhood, provided there is no 
other C at the new node yet. The rule reflects the fact that in the real 
life there is usually a close correspondence between food reserves of an 
individual and its breeding success. 

If a P appears at a node occupied by a C, an analogous situation takes 
place. P feeds on C and then produces an offspring with a probability 
bp. 

If there is no i? at a node occupied by a C, then the consumer dies 
with a probabihty dc. Its body D stays at the node until used by an 
R. 

Similarly, if there is no C at a node populated by a P, the predator 
dies with a probability dp and turns into a detritus portion D. 
If an R meets a D at a node, it feeds and then produces offspring at 
empty nodes in its Moore neighborhood. The eaten detritus is removed 
from the lattice. While the upper level species C and P are allowed to 
produce at most one progeny after each breeding success, for the basal 
species R we assume that it may populate all empty sites (i.e. 8 nodes 
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State in time t 



State in time t + 1 



E — > E 

R — > Rot E 

C ^ C oi D 

RC ^ C i+C) 

P — y P oi D 

D — > D 

RP — > RP or RD or P 

CP ^ P (+P) 

RCP — > P l+P) 

RD — > R (+ max 8R) 

CD — > CD OT D 

RCD — y CD 

PD — > PD oi D 

RPD — > RP (+ max 8R) 

CPD — > PD (+F) 

RCPD — > RP (+P,+ max 8R) 



Table 1: Look up table in our model. E, R, C, P and D stand for empty node, resource, 
consumer, predator and dead individual (consumer or predator), respectively. The values 
in parentheses represent a possible outcome of an action in the neighborhood of a given 
site. 



at most) in the neighborhood after breading, each of which with the 
probabihty h^. This assumption comphes with the fact that the lower 



the trophic level the higher its productivity [27 . 
7. If there is no D at a node occupied by an i?, it dies with a probability 
cLr and is removed from the lattice. Note that there is no conversion 
into a D in this case. 

As it follows from the above rules, a single node of the lattice may be 
in one of 16 states. Depending of its actual state different actions may 
be performed according to the rules. Possible states and outcomes of those 
actions are summarized in Tabled] Looking at this table one can immediately 
see that the system reveals a very rich dynamics. That is why we made some 
further assumptions. First of all, we take only occurrences of D into account, 
not its actual quantity. As a consequence, if a consumer or a predator dies 
at a node already containing a D, then its content will remain one D. In 
other words, we simply assume that resources are able to absorb all detritus 
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they find. Thus, its quantity does not really matter. A further presumption 
is that the dead resources do not contribute to the nutrient pool. 

A living species in our model does not move on the lattice, it remains 
localized. The only way of invading new lattice sites is via proliferation. In 
the simplest version of the model the directions of proliferation are chosen at 
random. However, it is possible to put some intelligence into the behavior of 
the species. 

The impact of different proliferation strategies was already discussed in 



Ref. j25|. In that work the following strategies has been taken into account: 



• Random - a mixture of the randomness and exclusion, i.e. the afore- 
mentioned basic rules. 

• CtR - consumers C put their progeny only on a site occupied by R, i.e. 
the direction of their invasion is driven by food availability. If there is 
no R in the neighborhood of C, no offspring is produced. Other species 
follow the Random strategy. 

• CaP - consumers C put their progeny on a site not occupied by P. 
Other species follow the Random strategy. 

• CtRaP - this is a mixture of the above strategies, i.e. consumers C 
put their progeny on a site occupied by i?, if its not occupied by P at 
the same time. Other species follow the Random strategy. 

• PtC - predators P put their progeny only on a site occupied by C. 
Other species follow the Random strategy. 

• CaPtC - predators P put their progeny only on a site occupied by C, 
whereas C avoid sites occupied by P. Resources R follow the Random 
strategy. 

If not stated otherwise, the results presented in the next sections of this 
paper have been obtained with the Random strategy. However we will refer 
to other strategies at many points for the sake of comparison. 

3. Results 

With the assumptions mentioned in the previous section, the model has 
seven important parameters: birth and death rates of each species {bi,di) 
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{i = R, C, P) and the linear size of the lattice L. To reduce the number of 
parameters we arbitrarily fix all death rates at 0.01. 

We set the linear size L of the lattice equal to 200 throughout the simu- 
lations. It was already shown in Ref. j24j that for system sizes larger than 
200 the viability of the food web depends only weakly on L. Thus, L = 200 
constitutes a reasonable compromise between the computational efforts and 
quality of the results. 

All results presented in this section were obtained in Monte Carlo simu- 
lations with the Moore neighborhood (i.e. 8 neighboring sites) on the square 
lattice, because it turned out to be a good compromise between the richness 
of the resulting dynamics of the system and the computational requirements. 
Among other choices we have tried the von Neumann neighborhood with only 
4 neighboring sites and extended Moore-like ones with 12 and 20 sites. While 
the von Neumann neighborhood led to an absorbing state (see Section |3TT] for 
more details) independently of the actual model parameters, the extended 
neighborhoods yielded results which were qualitatively the same compared 
to the less computational demanding Moore case. 

Most simulations were performed up to 10^ Monte Carlo steps (MCS). 
This particular value turned out to be a reasonable choice as well, because 
changing it to higher values did not significantly modified the results. In 
other words, it is highly unlikely in our model that a system, which is alive 
at T = 10^, will die out afterwards. 

If required (e.g. for phase diagrams in Sec. 13. 4p . we averaged the results 
over 100 independent runs in order to get a good statistics. 

3.1. Asymptotic states 

Since the asymptotic states of the system have been already discussed in 



Refs. [2J] and |25[, we give here only a brief summary of the most important 
findings. 

In case of proliferation strategies Random, CtR and PtC (see Sec. [2] 
for definitions) the system ends up in one of two distinct asymptotic states 
depending on the particular values of the model parameters: in an absorbing 
one with all species being extinct or in a coexisting one, in which all species 
survive till the end of a simulation. 

Theoretically, one could expect that a state with species R and C having 
non-zero concentrations and P being extinct is possible as well. Such a 
state was not observed in the simulations. Thus the predators seem to be 
indispensable for the survival of the system. 
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Three other strategies considered here, i.e. CtRaP, CaP and CaPtC, end 
up in an absorbing state independently of the actual model parameters. The 
reason is as follows. The presence of the predators in the system is required 
to control the number of consumers. Otherwise, consumer population will 
grow rapidly and kill all available resources driving the entire system to 
extinction. According to Table [H only nodes in the states CP, RCP, CPD 
and RCPD can lead to the migration of the predators, which is needed 
to maintain their population. In general, there are three sources of those 
four states: (1) states resulting from initial conditions, (2) predators which 
put their progeny on nodes occupied by C's and (3) consumers which put 
their offspring on nodes occupied by P's. In the strategies that end up in 
a coexisting state, all those sources are present. However, in the CtRaP, 
CaP and CaPtC strategies consumers avoid nodes occupied by predators. 



i.e. the third source is missing [25[. Although predators from initial CP 
and RCP states (source 1) are able to proliferate or even migrate farther if 
there are some consumers in their neighborhoods (source 2), there is no way 
for them to proceed with proliferation after all consumers in their vicinity 
will be eaten up because no P node will be populated by a new C. As a 
consequence, they will die out after a while. 

Closer look at the absorbing states reveals that there are two different 
ways for the system to reach them. Depending on the parameters of the 
model, in some cases the concentrations of species decrease almost monoton- 
ically towards zero. Other possibility is that well pronounced peaks in the 
i?, C and D densities appear before extinction. 

Examples of each asymptotic state in case of Random strategy are shown 
in Fig. |2j The curves were obtained from single runs. The parameters of the 
simulations may be found in the plots and their captions. 

We will show that the peaks seen in the bottom right plot in Fig. [2] 
correspond to density waves traveling through the system. Analysis of those 
waves is the main goal of this paper. 

3. 2. Origin of density outbursts 

In Fig. [3] the time evolution of R and C in the system from the bottom 
right plot of Fig [2] is presented again in a smaller time window. As far as the 
resources R are concerned, we observe rapid explosions of their population 
followed by almost equally rapid declines. This behavior bears a resemblance 



to characteristics of an excitable media 28 . Moreover, a similar time evolu- 



tion may be found in processes that take place in plankton populations and 
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bp = 0.8, bQ ^ 0.8, bp = 0.6 




1000 2000 3000 4000 5000 500 1000 1500 2000 2500 

Time (MCS) Time (l«CS) 



Figure 2: (Color online) Example of asymptotic states in case of the Random strategy. In 
the alive state (top plot) all species survive till the end of the simulation. In the absorbing 
state (bottom plots) the entire system dies out. The are two distinct ways for the system 
to reach the absorbing state: almost monotonic decrease of the species densities (bottom 
left) and well pronounced peaks of the densities before the extinction (bottom right). Birth 
rates used in the simulations are indicated in the plots. Death rates were all set to 0.01 
and the linear size of the lattice was equal to 200. 
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bp - 0.8, bp - 0.8, bp - 0.4 



bp - 0.8, be = 0.8, bp - 0.4 





500 600 700 800 900 1000 1100 1200 1300 
Time {MCS) 



500 600 700 800 900 1000 1100 1200 1300 
Time (MCS) 



Figure 3: Well pronounced peaks in the densities of R (left plot) and C (right plot). See 
bottom right plot in Fig. [2] for more details. 



are known as "spring blooms" and "red tide" phenomena 29 



Let us now have a closer look at the system at a microscopic level. Snap- 
shots of the system at different MC steps are shown in Fig. HI At t = 1200 
MCS, i.e. just before the last peak in the density of R shown in Fig. [3] occurs, 
there is almost no activity in the system. Most of the nodes are occupied 
by a C or a D and there are only a few resources and predators. However, 
there is a small cluster of R in the top center part of the lattice. This cluster 
grows, splits into two parts and moves away from its origin, as it may be 
seen in the snapshot at t = 1210 MCS. Note, that it is followed by an even 
bigger cluster of C's. Both clusters continue to grow and spread over the 
entire system. This fact is reflected by the peaks in the densities as may be 
seen in Fig. [31 Finally the resources hit the boundaries of the system (or the 
front of another wave) and have no place to further escape from C's. They 
will be diminished very quickly by the consumer's wave. Then the consumers 
will partially die out due to the lack of food - the densities of both species 
decline rapidly. 

Since the traveling waves presented above were observed only for some 
particular values of the simulation parameters, i.e. the birth rates of the 
species at different trophic levels, it is to expect that special conditions are 
necessary for an R wave to emerge (see Fig. [5]). First, there must be an 
RP node, so that R is protected from being eaten by a C for some time 
interval. If P dies due to the lack of food, it provides nutrient for R, which 
then proliferates to the neighboring sites. Second, if the origin RP node was 
surrounded by D's, the progenies are able to feed immediately and produce 
offspring by themselves. In this way the initial cluster seen in the top left 
snapshot in Fig. H] is created. Third, the overall concentration of D must be 
high in order to feed the wave through the entire process. 
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Figure 4: (Color online) Snapshots of the system at different MC steps. At t — 1200 MCS, 
most of the nodes are in a C (yellow) or D (grey) state. There is only a small cluster 
of R (red) in the top center part of the lattice. This cluster grows, splits into two parts 
and moves away from its origin, as it may be seen in the snapshot at t ^ 1210 MCS. It 
is followed by even bigger cluster of C"s. Both clusters continue to grow and spread over 
the entire system. Finally the resources hit the boundaries of the system (or the front 
of another wave) and have no place to further escape from C's. They will be diminished 
very quickly by the consumer's wave. Then the consumers will partially die out due to 
the lack of food and again, most of the nodes will be in a C or D state. Parameters of the 
simulation: L = 200,6/? = 0.8,6c = 0.8,6p = 0.4. 
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Figure 5: Origin of R blooms. A resource R on an RP node is protected by a predator P 
from being eaten by a consumer. If P dies, it provides food for R and lets it proliferate. 
If the origin node is surrounded by the detritus, newborns get enough food to reproduce 
and an avalanche starts. 
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Figure 6: Origin of C blooms. If an R wave hits a node populated by a C, it gets food 
and starts to produce offspring. Since there is plenty of food for C newborns brought by 
the wave, they start to invade the lattice and a C wave emerges. 

As may be seen from Fig. HI the R wave is followed by a C bloom. An 
explanation for this phenomenon may be found in Fig. [HI The C avalanche is 
triggered by the front of an R wave hitting a node occupied by a consumer. 
Since in all proliferation strategies under consideration the species R follows 
a random behavior pattern, in the situation presented in Fig. [6] it will put 
an offspring on the C node as well. As a result, the consumer gets food and 
starts to proliferate. Since there is plenty of food brought by the R wave, 
the new consumers start to invade the lattice - a C wave emerges. 

3.3. Spanning dusters 

In the previous section it has been already mentioned that an RP node 
and detritus states in its vicinity are necessary to feed the wave front (see 
Fig. E]). Let us now investigate these conditions in more detail. 
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Figure 7: Frequency of RP — > RD events overimposed on the resource density. Prolifera- 
tion strategy: Random. Parameters of the simulation: 6_r = 0.8, be = 0.8, bp = 0.4. 



We begin with checking the frequency of RP — > RD events. As it follows 
from Fig. [TJ although their frequency decreases slightly in the course of time, 
such events occur in the low density phase not only just before the emergence 
of a new wave. Thus, an RP node surrounded by the detritus is not a 
sufficient condition for the wave. 

In Fig. [HI distributions of detritus clusters at two different MC steps are 
presented. Both distribution were measured in the low density state of the 
resources. However, while the left plot falls into the time interval just after 
the density decline, the right plot describes a phase before the formation of 
a new wave. We see that the two distributions differ significantly. At the 
beginning of the low density phase, there are plenty of small detritus clusters. 
If an RP node is surrounded by such a small cluster, the emerging wave is 
suppressed sooner or later due to the lack of food for the R newborns. The 
situation changes upon entering the state corresponding to the right plot 
in Fig. [HI There are just a few small clusters and a big one with the size 
comparable with the size of the lattice. If an R wave starts to form at the 
boundary or in the bulk of the big cluster, it can then travel through the 
entire system. Of course, the cluster itself will be destroyed by the wave and 
it will take a while to rebuild it. That is why not all RP RD events lead 
to a density bloom and why we observe low density states in the system. 

Our hitherto findings indicate that there could be a connection between 
the emergence of density waves and the percolation of detritus on the lattice. 
In order to verify this hypothesis we checked in every step of our simulations 
if there exists a cluster of detritus spanning the whole lattice. The results are 
shown in Fig. [9l We see that indeed there are percolation periods (indicated 
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t = 1 1 40 [MCS] t = 1 200 [MCSJ 



Figure 8: Top: Detritus cluster distributions at two different MC steps. Parameters of the 
simulation: bji = 0.8, be = 0.8, bp = 0.4. At t = 1140 MCS (left plot), there are many 
small clusters on the lattice. The cluster distribution at t = 1200 MCS is different - there 
is one big cluster comparable with the size of the lattice and only few small ones. Bottom: 
Snapshots of D clusters corresponding to the above distributions. 
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600 700 800 900 1000 1100 
Time [MCS] 

Figure 9; Time evolution of R vs percolation intervals of D. To plot the intervals, in our 
simulations we assigned the value of 10 to each time step in which a detritus spanning 
cluster was detected, and otherwise. Parameters of the simulation: bji ~ 0.8, be = 0.8, 
bp = 0.4. 



by the rectangles) separated by time intervals, in which no spanning clusters 
occur. Moreover, the peaks in the resource density accrue at the ends of 
the percolation intervals. The reason is simply that after the emergence of a 
spanning cluster it takes some time for the wave to form and then to travel. 
It should be mentioned that no spanning clusters were observed in case of a 
coexistence state or an almost monotonic decline (see Section EH] for details). 
In both cases in every time step the cluster distributions were similar to the 
left plot of Fig. IHl It seems that the density blooms are possible only if the 
percolation of the detritus occurs in the system. 

Percolating detritus clusters shown in Fig. |9] may be interpreted as an ac- 
cumulation of nutrients in the system. In ecology, a response of an ecosystem 



to increased levels of nutrients is known as eutrophication [30|. High level 
of nutrients stimulates the primary production, causing a quick population 
increase (i.e. a bloom) of species such as algae in aquatic systems. These 
blooms may have many ecological effects, among others a decreased biodi- 



versity [3ll|. changes in species composition and dominance j32| and toxicity 
effects |33[. 

3.4- Phase diagrams 

To collect more information on the behavior of the system we performed 
a series of simulations with different sets of the model parameters. All birth 
rates of the species were varied from 0.1 to 0.9 with step 0.1. For each 
parameter set we run 100 independent runs. A set was tagged as an 'alive' one 
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Figure 10: A discrete 3D phase diagram in {bji, be, bp) space obtained for Random strat- 
egy. Circles correspond to the wave states, crosses are the ahve states. All birth rates take 
values from 0.1 to 0.9, with the step 0.1. 



if at least one out of those 100 runs ended up in a coexisting state. Similarly, 
a set was labeled as a 'wave' one if at least in one of the runs the density 
waves were observed. The collection of all 'alive' sets constitutes the 'alive' 
phase. The remaining sets form the 'dead' phase, which is divided into two 
subphases: the (almost) 'monotonic' one and the 'wave' one (see Section [XT] 
for more details). From our experience it follows that the transitions between 
phases are very sharp, i.e. if for a particular set of parameters one run ended 
in a coexistence state, then the most of the 100 runs led to the coexistence 
state as well. 

From the results of the above simulations one could generate a discrete 
3D phase diagram in {bji, be, bp) space. An example of such diagram is shown 
in Fig. HU] for the Random strategy. However, since those diagrams are not 
really readable, we will analyse its 2D sections at fixed resource's birth rates 
instead. 

In Fig. [m 2D phase diagrams in {be, bp) planes at br = 0.6 for different 
proliferation strategies are shown. For other values of bji we would get similar 
diagrams with only small changes in the areas of different phases. It is due to 



the fact that the viability of the system depends only weakly on bp \24l |25 



We see that there is a critical value of be for each proliferation strategy, 
above which the density waves appear in the system. In other words, only if 
the birth rate of the consumers is high enough, they are able to follow the 
resource's wave front according to the mechanism shown in Fig. El Otherwise 
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Figure 11: Phase diagrams in {be, bp) plane at bp — 0.6 for different proliferation strate- 
gies. Birth rates take values from 0.1 to 0.9, with the step 0.1. 
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Figure 12: Both total number of density wave states and total number of alive states in 
the 3D phase diagrams as functions of proliferation strategy. 



they will invade the lattice too slowly and the resources will die out before 
being caught up by the consumers. In this case we will observe after an initial 
density bloom the almost monotonic decline of the species densities towards 
the absorbing state (bottom left plot of Fig. [2]). 

Moreover, in the case of Random, CtR and PtR strategies, i.e. the strate- 
gies that support the 'alive' state as well, we also observe a critical value of 
bp, below which the density blooms occur. Thus the emergence of the waves 
require the birth rate of predators to be small enough to allow an almost un- 
controlled bloom of the consumers. If 6p is higher than the critical value, the 
predators are able to suppress the quick increase of the consumers' density 
and the system arrives at a coexistence state. These results corroborate the 
existence of "top-down" interactions that impart stability to food webs [14 . 

In order to provide more information on different proliferation strategies, 
we have counted both the 'alive' and 'wave' states in the discrete 3D phase 
diagrams as the one shown in Fig. [TOl The results are shown in Fig. [T2J 
Again, we see that the proliferation strategies split into two groups - those 
which support only the two absorbing states {CtRaP, CaP and CaPtC) 
and those which end up in a coexisting state as well {Random, CtR and 
PtC). As far as the first group is concerned, the strategies are very similar 
to each other, since the volume of the 'wave' phase depends only little on 
the particular strategy. Among the strategies building up the other group, 
there is a strong similarity between Random and CtR [25] . While the volume 
of the 'alive' phase is slightly bigger for the CtR strategy, we observe more 
'wave' states in the phase diagram of the Random one. The PtC strategy 
differs from the other two in the size of the phases. However, the differences 
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are smaller in case of the wave states. 



4. Conclusions 

A simple spatial model of a three level food web with a closed nutrient 
cycle has been investigated. The results complement our previous findings 



on the stability of such a system |24l . |25 



The time evolution of the model food web reveals two asymptotic states 
(Fig. 12]): an absorbing one with all species being extinct, and a coexisting one, 
in which concentrations of all species are non-zero. Theoretically, one could 
expect that a state with species R and C having non-zero concentrations 
and P being extinct is possible as well. Such a state was not observed in the 
simulations, because the predators are indispensable for the survival of the 
system. 

We found two possible ways for the system to reach the absorbing state 
depending on the particular values of the model parameters. In some cases 
the densities increase very quickly at the beginning of a simulation and then 
decline slowly and almost monotonically (small fluctuations disregarded). In 
others, well pronounced peaks in the R, C and D densities appear regularly 
before the extinction. 

Those peaks correspond to density waves traveling through the system, 
a phenomenon which is often observed in many complex systems in nature. 
Understanding the mechanism that triggers those waves was the focus of the 
present work. We have shown that several conditions have to be met for the 
waves to emerge: 

1. An RP node surrounded by the detritus is needed to initiate the re- 
source wave. 

2. A spanning cluster of the detritus on the lattice is required to feed the 
front of the resource wave. 

3. The consumers must be able to reproduce quickly to follow the resource 
wave. In other words, their birth rate be must be higher than a critical 
value. 

4. The birth rate bp of the predators must be lower than a critical value to 
keep the number of the predators low and to allow almost uncontrolled 
grow of consumers followed by their quick decline due to natural causes. 

The birth rate bji of the resources and the different proliferation strategies 
turned out to be less important for the emergence of the density waves. 
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The existence of a critical value of bp, above which the food web enters 
the coexistence phase (see Fig. [TT]) is consistent with the hypothesis of top- 
down interactions being essential for the stability of food webs IJ]. Indeed, 
the emergence of the waves require the birth rate of predators to be small 
enough to allow an almost uncontrolled bloom of the consumers. If bp is 
higher than the critical value, the predators are able to suppress the quick 
increase of the consumers' density and the system survives. 

The coupling between the percolation and spatio-temporal patterns has 



been already assumed or observed in many complex systems [34 



35, 36, 37 



Since the percolation of the detritus in our model corresponds to an accu- 
mulation of nutrients, the processes we observe in the model reveal some 



resemblance to the eutrophication phenomena [30| occuring in both aquatic 
and terrestrial ecosystems. Within our simple model this accumulation in- 
duces density blooms which drive the entire food web to extinction. In real 
ecosystems the outcome of eutrophication is usually not lethal for a system 
as a whole. Nevertheless, its ecological impact is severe, since it may result 
in loss of species, decreased biodiversity, changes in species composition and 
dominance or toxicity effects. 

Identifying percolation as one of the triggers for the density blooms may 
have a practical impact for the biological control. From the theory of perco- 



lation of disordered media [38| it follows that it is enough to destroy some 
critical fraction of "transmission-promoting" sites to suppress the propaga- 
tion of waves and their negative effects on the system. However, it should be 
checked which asymptotic state the system reaches after such an intervention. 
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